#Long formatted R1 data for all scans + frontal lobe regions at each depth
myelin.compartments.7T.long <- lapply(myelin.compartments.7T, function(depth){
  cols_to_pivot <- names(depth)[grepl("ROI", names(depth))]
  depth.long <- depth %>% pivot_longer(cols = cols_to_pivot, names_to = "orig_parcelname", values_to = "R1")
  depth.long <- merge(depth.long, glasser.frontal, by = "orig_parcelname")
  depth.long <- depth.long[!(depth.long$orig_parcelname %in% glasser.snr.exclude$orig_parcelname),]
  depth.long <- depth.long %>% mutate(subses = sprintf("%s_%s", subject_id, session_id))
  return(depth.long)
})

#Calculate average frontal lobe R1 for each sub/ses
myelin.frontallobe.7T <- lapply(myelin.compartments.7T.long, function(depth){
    depth <- depth %>% group_by(subses) %>% do(mean.R1 = mean(.$R1)) %>% unnest(cols = c( "mean.R1"))
    depth <- merge(depth, participants, by = "subses")
    depth$subject_id <- as.factor(depth$subject_id)
    depth$sex <- as.factor(depth$sex)
    return(depth)
})

#Function to fit a frontal lobe gam for a given depth (superficial, middle, deep) and return model summary
gam.summary <- function(depth){
  
  depth.data <- depth
  
  model <- gam(mean.R1 ~ s(age, k = 4, fx = F) + s(subject_id, bs = "re"), method = c("REML"), data = depth.data)
  
  print(summary(model))
}

#Function to fit a frontal lobe gam for a given depth and plot the developmental trajectory
plot.trajectory <- function(depth, display_color, y_lims, y_ticks){
  
  depth.data <- depth
  
  #First fit the gam and get fitted values and derivatives 
  depth.model <- gam.statistics.smooths(input.df = depth.data, region = "mean.R1", smooth_var = "age", id_var = "subject_id", covariates = "NA", random_intercepts = T, random_slopes = F, knots = 4, set_fx = F)

  #Age spline plot with participant-level data
  trajectory.plot <- ggplot() +
    geom_point(data = depth.data, aes(x = age, y = mean.R1, group = subject_id), color = "black", size = .4) +
    geom_line(data = depth.data, aes(x = age, y = mean.R1, group = subject_id), linewidth = 0.3, color = "gray75") +
    geom_line(data = depth.model$gam.fittedvalues, aes(x = age, y = fitted), color = display_color, linewidth = 1.5) +
    labs(x="\nAge", y=sprintf("R1\n")) +
    theme_classic() +
    scale_y_continuous(limits = y_lims, breaks = y_ticks) + 
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))
  
  return(trajectory.plot)
}

#Function to fit a frontal lobe gam for a given depth and return the fitted values df
depth.fittedvalues <- function(depth){
  
  depth.data <- depth
  
  #First fit the gam and get fitted values
  depth.model <- gam.statistics.smooths(input.df = depth.data, region = "mean.R1", smooth_var = "age", id_var = "subject_id", covariates = "NA", random_intercepts = T, random_slopes = F, knots = 4, set_fx = F)
  
  depth.fitted <- depth.model$gam.fittedvalues
  return(depth.fitted)
}

#Compute longitudinal change in frontal lobe R1 between paired longitudinal sessions, at each cortical depth
compute_R1_longitudinalchange <- function(input.df){
  
  input.df$mp2rage.date <- ymd(input.df$mp2rage.date)
  
  #session-to-session raw longitudinal change in R1
  R1.longitudinal.change <- input.df %>% group_by(subject_id) %>%
    mutate(R1.change.raw = mean.R1 - lag(mean.R1)) %>% 
    ungroup() %>% select(subject_id, session_id, age, mp2rage.session_number, mean.R1, R1.change.raw)
  
  #session-to-session longitudinal change in age
  R1.longitudinal.change <- R1.longitudinal.change %>% group_by(subject_id) %>%
    mutate(age.change = age - lag(age)) %>% 
    ungroup() %>% select(subject_id, session_id, age, mp2rage.session_number, mean.R1, R1.change.raw, age.change)
  
  #midpoint between-session age
  R1.longitudinal.change <- R1.longitudinal.change %>% group_by(subject_id) %>%
    mutate(mean.age = (age + lag(age))/2) %>% 
    ungroup() %>% select(subject_id, session_id, age, mp2rage.session_number, mean.R1, R1.change.raw, age.change, mean.age)
  
  #longitudinal change in R1 scaled by change in age (delta R1/ delta age = R1 change slope!)
  R1.longitudinal.change <- R1.longitudinal.change %>% mutate(R1.change.slope = R1.change.raw/age.change)
  
  #remove subjects with a single timepoint only
  R1.longitudinal.change <- R1.longitudinal.change %>% group_by(subject_id) %>% filter(n() > 1) %>% ungroup() 
}

myelin.frontallobe.longchange <- lapply(myelin.frontallobe.7T, function(depth){
  compute_R1_longitudinalchange(input.df = depth)
})

# >80% of longitudinal visits show an increase in frontal cortex R1 in superficial, middle, and deep cortex
lapply(myelin.frontallobe.longchange, function(depth){
  R1.longincrease.percent <- round((sum(depth$R1.change.raw > 0, na.rm = T) / sum(!(is.na(depth$R1.change.raw))))*100, 2)
  sprintf("%s percent of longitudinal visits show an increase in R1", R1.longincrease.percent)
})

#Plot within-person R1 + slope at younger and older scans 
myelin.frontallobe.longchange.alldepths <- do.call(rbind, myelin.frontallobe.longchange) #long df with all depths
myelin.frontallobe.longchange.alldepths <- myelin.frontallobe.longchange.alldepths %>% mutate(depth = factor(str_remove(row.names(myelin.frontallobe.longchange.alldepths), "\\..*"), ordered = T, levels = ordered_depths))

scans.t1t2 <- myelin.frontallobe.longchange.alldepths %>% filter(mp2rage.session_number < 3) #first get paired visits 1 and 2
scans.t1t2$mp2rage.session_number[scans.t1t2$mp2rage.session_number == 1] <- "Younger" 
scans.t1t2$mp2rage.session_number[scans.t1t2$mp2rage.session_number == 2] <- "Older"
scans.t1t2$mp2rage.session_number <- factor(scans.t1t2$mp2rage.session_number, levels = c("Younger", "Older"))

scans.t2t3 <- myelin.frontallobe.longchange.alldepths %>% filter(mp2rage.session_number > 1) #next get paired visits 2 and 3
scans.t2t3 <- scans.t2t3 %>% group_by(subject_id, depth) %>% filter(n() > 1) %>%  ungroup() # make sure we only plot paired data
scans.t2t3$mp2rage.session_number[scans.t2t3$mp2rage.session_number == 2] <- "Younger"
scans.t2t3$mp2rage.session_number[scans.t2t3$mp2rage.session_number == 3] <- "Older"
scans.t2t3$mp2rage.session_number <- factor(scans.t2t3$mp2rage.session_number, levels = c("Younger", "Older"))


ggplot() +
  geom_point(data = scans.t1t2, aes(x = mp2rage.session_number, y = mean.R1, color = depth)) +
  geom_line(data = scans.t1t2, aes(x = mp2rage.session_number, y = mean.R1, group = interaction(subject_id, depth), color = depth), linewidth = 0.25) + 
  geom_point(data = scans.t2t3, aes(x = mp2rage.session_number, y = mean.R1, color = depth)) +
  geom_line(data = scans.t2t3, aes(x = mp2rage.session_number, y = mean.R1, group = interaction(subject_id, depth), color = depth), linewidth = 0.25) + 
  theme_classic() +
  scale_x_discrete(expand = c(0.03, 0.03)) +
  scale_color_manual(values = c(alpha(depth_colorbar[1], 0.7), alpha(depth_colorbar[2], 0.7), alpha(depth_colorbar[3], 0.7))) +
  xlab("\n") +
  ylab("Frontal Cortex R1\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure3_supplementaryb/Longitudinal_slopes.pdf", device = "pdf", dpi = 500, width = 2.65, height = 3.2)

myelin.frontallobe.longchange.alldepths %>% 
  ggplot(aes(x = R1.change.slope, y = depth, fill = depth)) +
  geom_violin(aes(fill = depth), color = "white", adjust = 1.25) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  xlab("\nR1 Rate of Change (slope)") +
  ylab("Cortical Compartment\n") +
  scale_x_continuous(limits = c(-0.008, 0.022)) +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure3_supplementaryb/R1slope_depthplot_violin.pdf", device = "pdf", dpi = 500, width = 2.6, height = 3.2)

R1age.slope.longitudinal <- lapply(myelin.frontallobe.longchange, function(depth){
  R1.age.slope <- mean(depth$R1.change.slope, na.rm = T)
  return(R1.age.slope)
})

R1age.slope.longitudinal <- do.call(rbind, R1age.slope.longitudinal) %>% as.data.frame() %>% set_names("R1.change.slope")
R1age.slope.longitudinal$depth <- factor(row.names(R1age.slope.longitudinal), ordered = T, levels = ordered_depths)

R1age.slope.longitudinal %>% 
  ggplot(aes(x = R1.change.slope, y = depth, fill = depth)) +
  geom_point(shape = 23, size = 5, color = "white") +
  theme_classic() +
  scale_fill_manual(values = depth_colorbar) +
  scale_x_continuous(limits = c(0.004, 0.0056)) +
  xlab("\nR1 Rate of Change (slope)") +
  ylab("Cortical Compartment\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure3_supplementaryb/R1slope_depthplot_meanpoint.pdf", device = "pdf", dpi = 500, width = 2.15, height = 1.7)

myelin.frontallobe.longchange.alldepths <- myelin.frontallobe.longchange.alldepths %>% filter(!is.na(R1.change.slope)) #75 longitudinal visits for superficial, middle, deep
myelin.frontallobe.longchange.alldepths$depth <- as.numeric(myelin.frontallobe.longchange.alldepths$depth)
myelin.frontallobe.longchange.alldepths <- as.data.frame(myelin.frontallobe.longchange.alldepths)

summary(gam(R1.change.slope ~ s(mean.age, k = 4, fx = F) + s(depth, bs = "re"), method = c("REML"), data = myelin.frontallobe.longchange.alldepths))

R1rate.by.age.model <- gam.statistics.smooths(input.df = myelin.frontallobe.longchange.alldepths, region = "R1.change.slope", smooth_var = "mean.age", id_var = "depth", covariates = "NA", random_intercepts = TRUE, random_slopes = FALSE, knots = 4, set_fx = FALSE)

ggplot(R1rate.by.age.model$gam.fittedvalues, aes(x = mean.age, y = fitted)) +
  geom_line(color = "black", linewidth = 1) + 
  theme_classic() +
  scale_y_continuous(limits = c(0.0039, 0.0064), breaks = c(0.004, 0.005, 0.006)) +
  xlab("Mean Age") + 
  ylab("R1 Rate of Change (predicted)\n") +
  theme(
  legend.position = "none",
  axis.text = element_text(size = 6, family = "Arial", color = c("black")),
  axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
  axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
  axis.line = element_line(linewidth = .2), 
  axis.ticks = element_line(linewidth = .2)) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure3_supplementaryb/R1slope_byage_trajectory.pdf", device = "pdf", dpi = 500, width = 2.07, height = 1.57)

fs.qc <- arrow::read_parquet("/Volumes/Hera/Projects/corticalmyelin_development/BIDS/derivatives/surface_metrics/7T_BrainMechanisms_brainmeasures.parquet") %>% select(subject_id, session_id, lh_euler, rh_euler) #Euler numbers from fstabulate
fs.qc$session_id <- gsub('.{15}$', '', fs.qc$session_id) 
fs.qc$euler_number <- fs.qc %>% select(lh_euler, rh_euler) %>% rowMeans() 
fs.qc <- left_join(participants, fs.qc)

myelin.frontallobe.longchange.alldepths <- merge(myelin.frontallobe.longchange.alldepths, fs.qc)

change.positive <- myelin.frontallobe.longchange.alldepths %>% filter(R1.change.slope > 0)
change.negative <- myelin.frontallobe.longchange.alldepths %>% filter(R1.change.slope < 0)
t.test(change.positive$euler_number, change.negative$euler_number)

t.test(change.positive$mean.age, change.negative$mean.age)

t.test(change.positive$age.change, change.negative$age.change)

table(change.positive$sex) #53% female

table(change.negative$sex) #76% female

sex_data <- table(
  group = c(rep("Group1", nrow(change.positive)), rep("Group2", nrow(change.negative))),
  sex   = c(change.positive$sex, change.negative$sex)
)
chisq.test(sex_data)

#Function to calculate the number of significant age effects post-FDR
ageeffect.results <- function(depth){
  ageeffects <- depth
  n_roi <- nrow(ageeffects)
 
  #significant smooth terms
  aageeffects <- ageeffects %>% mutate(significant = p.adjust(ageeffects$GAM.smooth.pvalue, method = c("fdr")) < 0.05)
  sig_age <- aageeffects %>% filter(significant == TRUE) %>% nrow()
  sig_age_percent <- round((sig_age/n_roi)*100, 2)
                         
  sprintf("%s percent of regions show a positive significant developmental change in R1", sig_age_percent)
}

lapply(gam.statistics.alldepths, function(depth){
  ageeffect.results(depth)
})

gam.derivatives.alldepths.long <- do.call(rbind, gam.derivatives.alldepths)
gam.derivatives.alldepths.long <- gam.derivatives.alldepths.long %>% mutate(depth = factor(str_remove(row.names(gam.derivatives.alldepths.long), "\\..*")))

gam.derivatives.alldepths.long <- merge(gam.derivatives.alldepths.long, glasser.regions, sort = F)

#Calculate the average first derivative in each region at each depth
regional.rate.alldepths <- gam.derivatives.alldepths.long %>% group_by(label, depth) %>% do(rate = mean(.$derivative)) %>% unnest(cols = "rate")

#Calculate the average across-region first derivative at each depth
frontallobe.rate.alldepths <- regional.rate.alldepths %>% group_by(depth) %>% do(mean.rate = mean(.$rate)) %>% unnest(cols = "mean.rate")
frontallobe.rate.alldepths$depth <- factor(frontallobe.rate.alldepths$depth, ordered = T, levels = ordered_depths)

gam.statistics.alldepths.long <- do.call(rbind, gam.statistics.alldepths)
gam.statistics.alldepths.long <- gam.statistics.alldepths.long %>% mutate(depth = factor(str_remove(row.names(gam.statistics.alldepths.long), "\\..*")))

regional.maturation.alldepths <- gam.statistics.alldepths.long
regional.maturation.alldepths <- merge(regional.maturation.alldepths, glasser.regions, sort = F)

#Calculate the average across-region age of maturation at each depth
frontallobe.maturation.alldepths <- regional.maturation.alldepths %>% group_by(depth) %>% do(mean.age = mean(.$smooth.increase.offset, na.rm = T)) %>% unnest(cols = "mean.age")
frontallobe.maturation.alldepths$depth <- factor(frontallobe.maturation.alldepths$depth, ordered = T, levels = ordered_depths)

SA.frontal <- S.A.axis %>% filter(.$orig_parcelname %in% glasser.frontal$orig_parcelname,) %>% filter(!(.$orig_parcelname %in% glasser.snr.exclude$orig_parcelname)) %>% filter(label != "lh_L_10pp")

ggplot() + 
  geom_brain(data = SA.frontal, atlas = glasser, mapping = aes(fill = SA.axis), colour=I("gray20"), size=I(.08)) + theme_classic() + theme(legend.position = "none") + 
  scale_fill_gradientn(colors = map.colorbar, trans = 'reverse', limits = c(360, 60), oob = squish, na.value = "white") + 
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray50", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75")) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure4/SAaxis_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.3, height = 2.2)

rate.SAaxis.bydepth <- data.frame(depth = factor(), SA.corr = double())
for(this.depth in c(unique(regional.rate.alldepths$depth))){
  depth.data <- regional.rate.alldepths %>% filter(depth == this.depth)
  depth.data <- merge(depth.data, S.A.axis)
  corr <- cor(depth.data$SA.axis, depth.data$rate, method = c("spearman"), use = "complete.obs")
  corr <- data.frame(depth = this.depth, SA.corr = corr)
  rate.SAaxis.bydepth <- rbind(rate.SAaxis.bydepth, corr)
}
rate.SAaxis.bydepth

rate.SAaxis.bydepth$depth <- factor(rate.SAaxis.bydepth$depth, ordered = T, levels = ordered_depths)
ggplot(rate.SAaxis.bydepth, aes(x = SA.corr, y = depth,  fill = depth)) +
  geom_bar(stat = "identity") +
  theme_classic() +
  xlab("\nDevelopmental correlation") +
  ylab("Cortical compartment\n") +
  scale_x_continuous(breaks = c(-0.5, -0.25, 0), limits = c(-0.52, 0)) +
  scale_fill_manual(values = depth_colorbar) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure4/SAaxis_rate_depthcorrelation.pdf", device = "pdf", dpi = 500, width = 2.6, height = 1.7)

regional.rate.alldepths %>% filter(depth == "superficial") %>% merge(., S.A.axis) %>% 
  ggplot(., aes(x = SA.axis, y = rate)) + geom_point(size = .5) +
  geom_smooth(method = "lm", linewidth = 1.5, color = depth_colorbar[3], fill = "gray78") +
  labs(x="\nSensorimotor-association axis", y=sprintf("Rate of increase (derivative)\n")) +
  scale_y_continuous(breaks = c(0, 0.001, 0.002)) +
  theme_classic() +
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure4/SAaxis_rate_superficial.pdf", device = "pdf", dpi = 500, width = 2.6, height = 1.7)

bigbrain.frontal <- bigbrain %>% filter(.$orig_parcelname %in% glasser.frontal$orig_parcelname,) %>% filter(!(.$orig_parcelname %in% glasser.snr.exclude$orig_parcelname)) %>% filter(label != "lh_L_10pp")

ggplot() + 
  geom_brain(data = bigbrain.frontal, atlas = glasser, mapping = aes(fill = cytoarchitecture.gradient), colour=I("gray20"), size=I(.08)) + theme_classic() + theme(legend.position = "none") + 
  scale_fill_gradientn(colors = map.colorbar, limits = c(-.3, .3), oob = squish, na.value = "white") + 
  theme(legend.position = "none", axis.text = element_blank(), axis.line = element_blank(), axis.ticks = element_blank()) +
  new_scale_fill() + 
  geom_brain(data = glasser.plotting, atlas = glasser, mapping = aes(fill = cortex),  colour=I(alpha("gray50", 0)), size=I(0)) +
  scale_fill_manual(values = c(alpha("white", 0), "gray75")) 

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure4/Cytoarchitecture_corticalmap.pdf", device = "pdf", dpi = 500, width = 4.3, height = 2.2)

rate.bigbrain.bydepth <- data.frame(depth = factor(), histology.corr = double())
for(this.depth in c(unique(regional.rate.alldepths$depth))){
  depth.data <- regional.rate.alldepths %>% filter(depth == this.depth)
  depth.data <- merge(depth.data, bigbrain)
  corr <- cor(depth.data$cytoarchitecture.gradient, depth.data$rate, method = c("spearman"), use = "complete.obs")
  corr <- data.frame(depth = this.depth, histology.corr = corr)
  rate.bigbrain.bydepth <- rbind(rate.bigbrain.bydepth, corr)
}
rate.bigbrain.bydepth

rate.bigbrain.bydepth$depth <- factor(rate.bigbrain.bydepth$depth, ordered = T, levels = ordered_depths)
ggplot(rate.bigbrain.bydepth, aes(x = histology.corr, y = depth,  fill = depth)) +
  geom_bar(stat = "identity") +
  theme_classic() +
  xlab("\nDevelopmental correlation") +
  ylab("Cortical compartment\n") +
  scale_x_continuous(breaks = c(0, 0.35, 0.7)) +
  scale_fill_manual(values = depth_colorbar) +
  theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure4/Cytoarchitecture_rate_depthcorrelation.pdf", device = "pdf", dpi = 500, width = 2.6, height = 1.7)

regional.rate.alldepths %>% filter(depth == "superficial") %>% merge(., bigbrain) %>% 
  ggplot(., aes(x = cytoarchitecture.gradient, y = rate)) + geom_point(size = .5) +
  geom_smooth(method = "lm", linewidth = 1.5, color = depth_colorbar[3], fill = "gray78") +
  labs(x="\nCytoarchitectural variation", y=sprintf("Rate of increase (derivative)\n")) +
  scale_y_continuous(breaks = c(0, 0.001, 0.002)) +
  theme_classic() +
    theme(
        legend.position = "none", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text = element_text(size = 6, family = "Arial", color = c("black")),
        axis.title.x = element_text(size = 7, family ="Arial", color = c("black")),
        axis.title.y = element_text(size = 7, family ="Arial", color = c("black")),
        axis.line = element_line(linewidth = .2), 
        axis.ticks = element_line(linewidth = .2))

ggsave(filename = "/Users/valeriesydnor/Documents/ResearchProjects/CorticalMyelin_PlasticityDevelopment/Publications/Publication_Images/Figure4/Cytoarchitecture_rate_superficial.pdf", device = "pdf", dpi = 500, width = 2.6, height = 1.7)

regional.rate.alldepths <- merge(regional.rate.alldepths, S.A.axis)
regional.rate.alldepths <- merge(regional.rate.alldepths, bigbrain)
regional.rate.superficial <- regional.rate.alldepths %>% filter(depth == "superficial")
regional.rate.middle <- regional.rate.alldepths %>% filter(depth == "middle")
regional.rate.deep <- regional.rate.alldepths %>% filter(depth == "deep")

summary(lm(rate ~ SA.axis + cytoarchitecture.gradient, data = regional.rate.superficial))

summary(lm(rate ~ SA.axis + cytoarchitecture.gradient, data = regional.rate.middle))

summary(lm(rate ~ SA.axis + cytoarchitecture.gradient, data = regional.rate.deep))

gam.smoothestimates.alldepths.long <- do.call(rbind, gam.smoothestimates.alldepths)
gam.smoothestimates.alldepths.long <- gam.smoothestimates.alldepths.long %>% mutate(depth = factor(str_remove(row.names(gam.smoothestimates.alldepths.long), "\\..*")))